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One-body Green’s function theories implemented on the real frequency axis offer a natural for¬ 
malism for the unbiased theoretical determination of quasiparticle spectra in molecules and solids. 
Self-consistent Green’s function methods employing the imaginary axis formalism on the other hand 
can benefit from the iterative implicit resummation of higher order diagrams that are not included 
when only the first iteration is performed. Unfortunately, the imaginary axis Green’s function does 
not give direct access to the desired quasiparticle spectra, which undermines its utility. To this end 
we investigate how reliably one can calculate quasiparticle spectra from the Extended Koopmans’ 

Theorem (EKT) applied to the imaginary time Green’s function in a second order approximation 
(GE2). We find that EKT in conjunction with GE2 yields IPs and EAs that systematically un¬ 
derestimate experimental and accurate coupled-cluster reference values for a variety of molecules 
and atoms. This establishes that the EKT allows one to utilize the computational advantages of 
an imaginary axis implementation, while still being able to acquire real axis spectral properties. 

Because the EKT requires negligible computational effort, and can be used with a Green’s function 
from any level of theory, we conclude that it is a potentially very useful tool for the systematic study 
of quasiparticle spectra in realistic systems. 

I. INTRODUCTION 

Finding practical and principled methods to numeri¬ 
cally solve the many-electron Schrodinger equation for re¬ 
alistic chemical systems is a substantial problem that has 
been attacked by the scientific community for more than 
five decades now. Though diverse in approach, the bulk 
of these methods can be classified as either being based 
on the single-particle density p(r) or the many-body 
wavefunction Kohn-Sham density functional theory 
(DFT)[I1[2] and the coupled cluster theory (CC)[3]-[5] are 
two remarkably successful examples of these former and 
latter categories. Despite their merits and widespread 
use, by now it’s become well established these classes of 
methods have their respective limitations that are un¬ 
likely to be overcome in the near future. For instance 
CC theory, though being the gold standard for weakly 
correlated systems, cannot successfully account for multi¬ 
reference effects present in strongly correlated molecules 
and solids, and its computational scaling seems to limit 
use for large systems. DFT in contrast is computation¬ 
ally affordable [6], but is plagued by a lack of systematic 
approaches to improve functionals [7], as well as the fact 
that many approximate functionals can give a frustrat- 
ingly non-uniform performance across different systems 
and properties. This problem is exemplified for metal ox¬ 
ides, where even within the same system different DFT 
approximations are required for different properties [8l|9]. 

Furthermore, as spectral properties are concerned, even 
with the exact functional the Kohn-Sham (KS) gap can¬ 
not provide a theoretically sound approximation to the 
fundamental gap pUHIS] . 

Approaches based on the self-consistent single-particle 
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Green’s function, G{uj), are an attractive alternative with 
strengths and weaknesses that are complementary to den¬ 
sity and wavefunction-based methods. Similar to DFT, 
Green’s function methods can be implemented in a black¬ 
box manner in the atomic orbital (AO) basis pAUTT] . 
and typically will feature lower order polynomial scaling 
than wavefunction-based methods such as CCSD (cou¬ 
pled cluster singles doubles). Additionally, Green’s func¬ 
tion methods offer a natural language for embedding ap¬ 
proaches since the Green’s function can be partitioned 
among subsystems (similar to the density). 

Because of the need to evaluate G{uJn) on a numeri¬ 
cal grid, Green’s function based approaches can be 
broadly grouped into two classes: 

Greens funetions on the imaginary (Matsubara grid) 
axis describe a grand canonical ensemble and can be 
used to calculate temperature dependent quantities such 
as the free energy, specific heat, etc. A Green’s func¬ 
tion of the imaginary axis is a smooth function, with 
only a single pole near zero that is not accessed for 
any finite temperature larger than zero. Consequently, 
the imaginary axis formulation is a natural choice for 
any self-consistent approach where the Green’s function 
is expressed as a functional of the self-energy, G[Il(co’)], 
and where at self-consistency infinite classes of diagrams 
can be included due to the implicit resummation. A 
result of this diagrammatic resummation is that these 
approaches [T5j [161 EHl IlS] give finite results in strongly 
correlated cases when other methods such as truncated 
CC or MP2 (second order Mpller-Plesset [20 ] ) would di¬ 
verge pathologically. 

Greens funetions on the real frequeney axis are more 
commonly employed in quantum chemistry for zero- 
temperature calculations. The single-particle Green’s 
function on the real axis has multiple poles which cor¬ 
respond to ionization potential (IP) and electron affinity 
(FA) peaks [2114291 in the photoelectron spectrum. The 
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two-particle real-axis Green’s function is capable of de¬ 
scribing optical/neutral excitations. Performing the cal¬ 
culation self-consistently (where G[S(co’)] is evaluated as 
a functional of S(co’)) is notoriously difficult on the real 
axis, because the series of poles in the real-axis Green’s 
function requires a non-uniform grid that can change be¬ 
tween iterations. 

The real and imaginary axis formalisms should be 
treated as complementary but requiring the development 
of different tools in order to evaluate accurate Green’s 
functions. The interesting question that arises is if one 
can use an approach that has “the best of both worlds”; 
that is, to calculate the Green’s function on the imagi¬ 
nary axis in order to take advantage of the efficient self- 
consistency, and to subsequently employ the imaginary 
axis solution to calculate spectra on the real axis. 

One of the routes to obtaining a real axis Green’s func¬ 
tion from the imaginary axis data is through the process 
of analytic continuation j3QH32] . This procedure is known 
to suffer from several problems, namely it cannot recover 
sharp spectral features, frequently is problematic in re¬ 
covering fundamental gap edges, and is very sensitive to 
the initial imaginary axis data. 

Here, we attempt to examine an alternative answer 
to this problem. We will investigate how accurately 
and reliably one can calculate quasiparticle spectra from 
the extended Koopman’s theorem (EKT)|3lS5] start¬ 
ing from self-consistent Green’s function many-body 
theory [161 HI H] in an imaginary axis implementation. 
The EKT is valuable because it allows one to obtain, 
in principle, both IPs and EAs from a single electronic 
structure calculation on the neutral system. Similar to 
population analysis [46], the actual EKT procedure itself 
takes place in a simple post-calculation analysis, utilizing 
the density matrix and other quantities obtained by the 
preceding correlated method. As a result the EKT can be 
implemented in a blackbox manner, requires only a neg¬ 
ligible fraction of time, and is not tied to any particular 
level of theory. Despite its simplicity, most of the efforts 
so far seem to have focused on its application for calcu¬ 
lating IPs While an 

interesting benchmark of the EKT for EAs has appeared 
recently [45], in that work Bozkaya actually obtained the 
EA indirectly by calculating the IP of the anion. In con¬ 
trast, in this work we will calculate IPs and EAs from 
the neutral system via the EKT. As far as we are aware, 
EAs typically have not been calculated in this manner. 

While the EKT has usually been formulated in terms 
of a generalized Eockian[40] that is evaluated using meth¬ 
ods such as configuration interaction (Cl) or MP2, in this 
work we approach EKT with the machinery of second or¬ 
der Green’s function theory (GE2) [TblfTS] . As shown by 
Dahlen, Stan, and van Leeuwen[T8l [19], via EKT it is 
possible to calculate both IPs and EAs from the imagi¬ 
nary time Green’s function of the neutral system alone. 
This is potentially very useful for studying the spectral 
properties of extended systems, because it would circum¬ 
vent the numerically ill defined step of analytic continu¬ 


ation to the real axis. While their initial results for IPs 
obtained with GE2 and GW were promising, here we in¬ 
tend to examine GE2’s performance for IPs and EAs for 
a wider group of atoms and small molecules. 


II. THEORETICAL CALCULATIONS OF 
IONIZATION POTENTIALS AND ELECTRON 
AFFINITIES 

To make this work self-contained, here we briefly re¬ 
view some of the different strategies that have been uti¬ 
lized to calculate IPs, EAs, and the fundamental gap 
Eg = IV — EA, using wavefunction theory, DET, many- 
body theory, and combinations thereof. Naively, the 
simplest strategy for obtaining these quantities would 
be by energy differences of the neutral and charged 
cation/anion systems. However this is fundamentally 
problematic for periodic boundary-conditions (PBC) cal¬ 
culations of materials in the solid-state, where it would 
imply the presence of an infinite amount of unbalanced 
charge in the crystal. Even for finite systems, careful 
early Cl studies j49l - [5T] found that obtaining accurate 
EAs from energy-differences was particularly challenging 
because the correlation energy of the anion could con¬ 
verge appreciably slower than that of the neutral system. 
Additionally, methods with significant amounts of self¬ 
interaction error (such as approximate DET) may not be 
able to bind some anions at all[52H54]. 

Eor these reasons, significant interest has been placed 
in obtaining IPs and EAs directly from a single cal¬ 
culation on the neutral system. An exemplar for this 
is offered by the equation of motion coupled cluster 
theory fEOM-CCi [554[59| for electron attachment and 
removal [6QU66) . EOM-CC is a generalization of the orig¬ 
inal coupled cluster theory to charged and neutral ex¬ 
citations, and thus largely inherits the advantages and 
disadvantages of the CC method: if a single determinant 
description is valid, and the calculation not prohibitively 
expensive, one can expect very accurate results. Closely 
related to IP/EA-EOM-CC is the coupled cluster Green’s 
function method of Nooijen and Sniiders[67]. that again 
will bear some of the advantages and disadvantages of 
the underlying CC theory. 

In the same spirit. Green’s function (also called elec¬ 
tron propagator) methods such as the n^^-order Alge¬ 
braic Diagrammatic Approximation fADCfnP [24l[68l[69] 
and various self-energy approximations [70] have found 
regular use for the accurate, direct determination of IPs 
and EAs in finite systems [2^ [25] (TThSO] . Typically the 
underlying scheme in these methods is the iterative di- 
agonalization of the self-energy on the real axis in the 
molecular orbital (MO) basis, 5](A^)c^ = A^c^, until one 
converges to a given IP, A^. As such, this is quite dis¬ 
tinct from the fully self-consistent Green’s function im¬ 
plementations in the local AO basis that have appeared 
recent Iv fTSUTT] and that we are considering presently in 
this work. 
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A theoretically not fully justified however computa¬ 
tionally cheap strategy for directly obtaining quasiparti¬ 
cle spectra would be to simply perform a DFT calculation 
on the neutral system, and then interpret the resulting 
KS eigenvalues as IPs and EAs. Unfortunately, as previ¬ 
ously mentioned there is no rigorous theoretical ground¬ 
ing for such a procedure [TOl [TTl |T3l |HTJ [82]: even if the 
exact functional were used, while the highest occupied 
KS orbital eigenvalue would be the negative of the ion¬ 
ization potential, the KS gap Eg^ would still differ from 
the exact fundamental gap Eg by an amount equalling 
the derivative discontinuity of the exchange-correlation 
(XC) potential. Eg = Eg^ + A^c- As shown recently 
with stretched hydrogen chains, when strong correlation 
is present the derivative-discontinuity A^c can actually 
become the dominant contribution to Eg m- If one 
disregards this principled objection, in practice the KS 
gap still underpredicts the fundamental gap significantly, 
which is partly a consequence of the fractional-charge er¬ 
ror in approximate density functionals [831187] . 

Because Hartree-Fock and KS DFT tend to have op¬ 
posing fractional-charge errors, these methods will typ¬ 
ically overestimate and underestimate the fundamental 
gap, respectively. For this reason there will almost al¬ 
ways be a system-dependent empirical hybrid functional 
that will describe the fundamental gap well[8l|9l|88] when 
calculated in the Generalized Kohn-Sham (GKS)[89] for¬ 
malism as band energy-differences. From this stand¬ 
point it’s understandable why the HSE functional [90] 
has enjoyed so much success for many insulators and 
semiconductors [9TI - [93] . though it’s been remarked the 
HSE gap will tend to match the optical gap better than 
the fundamental gap[92l|93]. 

Because of these fundamental issues with the KS eigen¬ 
values, a very commonly used strategy has been to per¬ 
form “one-shot” corrections to the DFT spectra with the 
Go Ho approximation [m jOU [95]. From a purely prag¬ 
matic viewpoint the resulting quasiparticle spectra can 
be much improved with respect to experiment. However 
from a principled standpoint this is not entirely satisfy¬ 
ing because it is highly dependent on the starting DFT 
solution which is functional dependent. Because of this, 
the Go Ho results can vary significantly depending on the 
combination of system and density functional used [96]- 
unsj. To highlight one such example, for hematite {a- 
Fe 203 ) the GoH^o correction can yield quasiparticle gaps 
ranging from 1.3, 4.0, to 4.5 eV depending on whether 
PBE |lQ6j . HSE[^, or PBEQ [TQ7] is used, respectively, 
which is contrasted with the experimentally determined 
gap of 2.6T0.4 eV |lQQj . Though controversy exists in 
whether self-consistency will overall improve or worsen 
results m EH EHl [TQ81 ItT 2] . clearly fully self-consistent 
Green’s function calculations are valuable because they 
are reference independent. We review one such self- 
consistent implementation next. 


III. GF2 THEORY AND IMPLEMENTATION 

The real axis single-particle Green’s function, G(cc;), 
determines the expectation value of all single-particle ob¬ 
servables, in addition to the spectral density of states, 
IPs, and EAs. Unfortunately, calculating G{uj) exactly 
for a large system is not any more feasible than calculat¬ 
ing the exact wavefunction T. Nonetheless we can cal¬ 
culate the Green’s function of a non-interacting system, 
Go(cc;), very easily, and then correct it for the missing 
many-body correlation effects in a systematic way via 
the Dyson equation. Given some Go(co’), the exact G{uj) 
can be obtained by expanding in terms of the proper 
self-energy, F(co’), and analytically summing to yield the 
Dyson equation 

G{lj) = Gq{uj) + Go(u;)F(cj)Go(cj) 
+Go(co’)F(cj)Go(cj)F(c(j)Go(cj) + • • • 

^ n 

The self-energy, F(u;), is a frequency-dependent single¬ 
particle potential that encompasses all of the exchange- 
correlation (XG) effects of the many-body system. Anal¬ 
ogous to the KS potential Vxc^ one could think of H(uj) 
as being the XC potential that connects the Green’s func¬ 
tion of the non-interacting system to the Green’s func¬ 
tion of the fully interacting system. However, it is im¬ 
portant to remember that X(u;) is dynamic, nonlocal, 
and orbital-dependent, while Vxc in approximate DFT is 
typically static, “semilocal”, and density-dependent. In 
principle the exact E(u;) can be expanded diagrammat- 
ically in Gfu;j |113l 1114] . In a practical implementation, 
one chooses a subset of diagrams that can be evaluated 
in a computationally tractable manner. The resulting 
self-energy can then be written as an approximate func¬ 
tional of the Green’s function, F[G(ci;)], yielding a self- 
consistent set of equations. In this work we investigate 
the second order approximation (GF2), which includes 
all diagrams to second order and is shown in Fig. 

To take advantage of the easy to converge self- 
consistency procedure on the imaginary axis, the Green’s 
function can be written in a non-orthogonal atomic- 
orbital (AO) basis as 

GH =[(m + w)S-F-S(w)]“^ , (2) 

where S and F are the overlap and Fock matrices, fi is 
the chemical potential, uj is an imaginary frequency, and 
5](c(;) is the aforementioned frequency-dependent self¬ 
energy within the GF2 approximation containing second 
order diagrams from Fig. We use a uniform grid of 
imaginary Matsubara frequencies = (2n + l)i7Tll3^ 
with a power mesh imaginary time grid [11 5] running on 
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FIG. 1. Feynman diagrams included in the GF2 self-energy. 
A black arrow line represents a Green’s function, while a red 
squiggle represents a two-electron integral. The first two di¬ 
agrams are the frequency-independent Hartree and exchange 
terms, and are included in the Fock matrix. The next two di¬ 
agrams are the frequency-dependent pair-bubble and second 
order exchange terms, and represent the second order corre¬ 
lation effects covered by XI (cj). 


(i.e. /i ~ 2 ^^homo T ^lumo)-) F — ^HF-) s^nd 5](cj) — 0) 
generated via output from the Dalton electronic structure 
Drogram |116j . At self-consistency G(ci;) will not depend 
on the starting-reference [17], though practically some ini¬ 
tial guesses might be better than others for converging 
rapidly. As a final note, it should be understood that the 
Green’s function depends on /i, and therefore P implicitly 
depends on /i as well. This means the chemical poten¬ 
tial will need to be adjusted from iteration-to-iteration 
to maintain the correct electron number. 

For purposes of comparison, we will also consider a 
non self-consistent Green’s function obtained from the 
first iteration of the Dyson equation, given by 


the interval 0 < r < /3, where p is the inverse temper¬ 
ature. We choose to build the Green’s function on the 
frequency axis because of the simplicity of Eq. in con¬ 
trast to the expression for G(r) which is more cumber¬ 
some and requires integrations over r points [T8|. Once 
built, the imaginary frequency Green’s function can be 
fast Fourier transformed (FFT) to the imaginary time 
domain, G(r). The correlated density matrix P then is 
evaluated as 


P = -2G(r=/3) . (3) 

Provided with P, the correlated Fock matrix is built by 


Fij = hij + 'Y^Pkiiyijik - -j^ikij) , (4) 

kl 

where hij and Yijik are one and two-electron integrals 
in the AO basis. Note that the (frequency-independent) 
first-order self-energy is already covered by the Hartree- 
Fock mean-field, Soo = FkiFii^ijik - Finally, 

the (frequency dependent) second order self energy can 
be built in the time-domain as 


Sii(r) 


- E Gkl{T)G„^n{T)Gpq{-T) 
klmnpq 


^^imqk {^^Ipnj ^nplj^ , 


(5) 


and then FFT to the frequency domain. It is simpler 
to build X) on the time axis, because in the r domain 
the self-energy factorizes into simple products of Green’s 
functions, whereas in the uj domain it requires integra¬ 
tions of Green’s functions over frequencies. Furnished 
with an updated F and X)(cj), we can return to Eq.j^and 
rebuild G(cc;). Taken altogether Eq.[^[^|^ and [^present 
a self-consistent procedure for solving the Dyson equation 
in a second order approximation to the self energy. To 
initiate the self-consistency an approximate zeroth order 
Green’s function, Go(cj), is necessary, which practically 
can be supplied by DET or Hartree-Eock (HE) calcula¬ 
tions. In this work we use an initial HE Green’s function 


Gi(w)=[(/i + w)S-Ffff-S[Go(a;)]] ' . (6) 

Here X)[Go(co’)] is simply the self-energy obtained when 
Go, which in our case is Ghf^ is inserted into Eq. 
and /i is set so that Go(ci;) has good particle number. 
Eor conciseness, we will refer to this simply as GoE2, in 
analogy to Go Wo- Since this Green’s function is not self- 
consistent, it will carry a starting reference dependence. 

Since the Green’s function obtained by self-consistent 
or non-self-consistent GE2 is expressed on an imagi¬ 
nary grid, we aim to employ the Extended Koopman’s 
Theorem (EKT)[33l [34l [36l [571 ESI SH SS] to obtain 
ionization potentials and electron affinities, which are 
real axis quantities and can be used to produce the 
spectral density of states A{uj) expressed as A(u;) = 
— ^Tr[ImG(cc;)S]. We give a brief discussion of the EKT 
theory and implementation next. 


A. Extended Koopmans’ Theorem 


Given a system with Hamiltonian H and N electron 
state I at) satisfying H\n) = E\n), by using second- 
quantized operators, a\N) = cl^n) = |a^+i), the 

energies of the anion, neutral, and cation states can be 
expressed as 


En+1 = (jv|aiJa'l'|iv) , 

En = {n\H\n) , (7) 

Em -1 = {N\a^ Ha\N) , 

respectively. It should be understoood these operators 
are expanded in a basis, d = Ci(j)i^ with the expansion 
coefficients q chosen so that the anion (cation) state re¬ 
mains normalized [39]. Provided with Eq.j^the ionization 
potential (/) and electron affinity {A) can be defined as 


/ = En-1 - Em = -(iv|a’l'[a,iJ]|iv) , 

A = Em- Em+1 = (iv|a[atF]|iv) . 
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A Lagrangian for / and A can now be constructed, given 
by 


Cl =-{N\a'‘[a,H]\N)+e{{N\a^a\N)-1) , 

Ca = (jv|a[a^ ff]|jv) + e((iv|aa^|jv) - 1) , 

where the right-hand term constrains the cation/anion 
state to be normalized. Expanding the operators in their 
basis, and exploiting that aa^ = 1 — d’^'d, the station¬ 
ary solution 5CI5c\ = 0 of Eq. [^yields the generalized 
eigenvalue problem 


H^c = ePc , 

= eP^c , 


( 10 ) 


where ^, H]\n) and = 

— {N\^i[^pH]\N) are generalized Eock matrices, and 

[P]ij = (AT I 0^1 at) is the density matrix. P^ is the vir¬ 
tual (or “hole”) density matrix, which is defined within 
the orthogonal Lowdin basis as P^ = 21 — P, where I 
is the identity matrix and the factor of two accounts for 
double occupation in our spin-restricted formalism, or 
equivalently in terms of Green’s functions in the Lowdin 
basis as P^ = — 2 G(r=o+). 

Eor practical calculations we want to connect Eq.p!Q]to 
Green’s function many-body theory in the following way: 
Erom the definition of the time-dependent single-particle 
Green’s function |114j . G(r), one can show 


dG^J{N\a[a^,H]\N) , r >0 
r^o dr ^(A^|d'*'[d, ^]|A^) , r <0 


( 11 ) 


Eor compactness we simply write these two possibili¬ 
ties as 9^G(r) |o+ and 9^G(r)|o-. Eq. has this form 
because of the discontinuity in the Green’s function at 
r = 0 , and furthermore G(r) = —G(r -h /d), which means 
that 9T-G(r)|o- = —drG{T)\j 3 . The important point of 
Eq. is the generalized Eockians appearing in the stan¬ 
dard EKT Eq. can be replaced with time-derivatives 
of the Green’s function on the imaginary-domain in the 
AO basis. Introducing a matrix representation for the 
Green’s function in this basis, G(r), performing the 
transformation c = P“^Gc' (or c = P^; and multi¬ 

plying on the left with P“^/^ (or P^ ^^^) results in 


our spin-restricted formalism. We emphasize Eq. as¬ 
sumes that G(r) and P have been pre-transformed to the 
Lowdin basis. Eq. shows that diagonalization of the 
A_ and A+ matrices gives eigenvalues which, after sub¬ 
tracting out the chemical potential /r, yield the ionization 
potentials and electron affinities respectively. Eor exam¬ 
ple, if the Hartree-Eock Green’s functions were inserted 
in Eq. then this would yield simply the Koopman’s 
theorem IPs and EAs (in fact this is a good way to check 
the accuracy of one’s grid). Conceptually, the form of 
can be understood by realizing that G(r) for r 
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Eq. 

near (3 describes the particle distribution of the system 
and consequently electron removal, whereas for r near 
0 + it contains information on the hole distribution and 
therefore electron attachment. 

One subtlety is that diagonalizing A± will of course 
yield as many eigenvalues as there are AO basis functions, 
A 5 , yet only some of these eigenvalues may be physically 
meaningful as IPs or EAs. We find the simplest way to 
identify the correct eigenvalues is by the corresponding 
Dyson occupations 


PC_ =D , 

t “ (13) 

C^P„C+ = D„ . 

Here C± is the matrix of eigenvectors, C = {c^^, C 2 ,.. 
obtained from diagonalizing A±. The diagonal ele¬ 
ments of D (D^) correspond to occupations of Dyson 
orbitals for electron removal (attachment). As a con¬ 
sistency check one should find that Tr[D] = A, and 
Tr[D^] = 2 A 5 — N. Essentially the orbitals with large 
occupations (roughly speaking I < [D]^^ < 2) will indi¬ 
cate the IPs/E As one is interested in. 

As a final note, we stress again that the use of the 
Extended Koopmans’ Theorem is not limited to GE 2 or 
even Green’s function methods, and has been employed 
with a variety of methods at different levels of theory in 
the past [IHl[191 ES|34l [371 HOI 1111 IlllIlH]. It has been a 
matter of debate whether or not the lowest IP given using 
EKT is exact, or whether or not higher IPs obtained from 
this method are physically meaningful [40l|43]. We do not 
investigate this in this paper, but rather, show that EKT 
offers reasonable values for IPs and EAs. Presumably, if 
one were to use a more accurate Green’s function, one 
would obtain more accurate IPs and EAs. 


IV. COMPUTATIONAL DETAILS 


A_c' = e/c' , 6/ = / + /i 
A+c' = e^c' , ca = A + /i 
A_ =2P-i/2a,G(r)|o-P-'/2 (12) 

A+ = -2P;i/25^G(r)|o+P„-'/2 

where P = — 2G(r=/3), and P^ = — 2G(r=o+). The fac¬ 
tor of two in Eq. accounts for double occupation in 


Our GE 2 and GiE 2 calculations were carried out on 
an imaginary grid with 20,000 frequency points and 
4,400 time points[TT5], with an inverse-temperature of 
13= 100.0 [a.u~^]. Experimentally determined geome¬ 
tries were used for the molecules |117j . The basis sets 
used were Dunning’s aug-cc-pVXZ series [TTSl I119j . Re¬ 
stricted Hartree-Eock calculations carried out in the 
Dalton program [116] were used to generate an initial 
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FIG. 2. Ionization potentials and electron affinities of atoms calculated with the aug-cc-pVXZ basis set. 


Green’s function Go(co’) as input for our GF2 proce¬ 
dure. All GF2 and GiF2 calculations reported here are 
all-electron. The Hartree Fock IPs and EAs were ob¬ 
tained from standard Koopmans’ Theorem as the neg¬ 
ative of the HOMO and LUMO eigenvalues, respec¬ 
tively. To obtain accurate reference theoretical values, 
IPs and EAs were computed from energy differences 
of the charged and neutral species, using all-electron 
unrestricted UCCSD(T) calculations with the Gaussian 


09 package |12Qj (“UCCSD(T)=Eun” keyword). Let us 
stress that since IP and EA values in UCCSD(T) are cal¬ 
culated as a difference between the charged and neutral 
species, the obtained values not only include the bene¬ 
fit of energy lowering due to the use of an unrestricted 
method but also can take advantage of error cancella¬ 
tions. This stands in stark contrast to the IPs and EAs 
calculated from GE2 that is based on a restricted refer¬ 
ence (RHE) and does not benefit from error cancellation 
due to calculating differences. 


V. RESULTS 


Our ability to obtain accurate IPs and EAs will be af¬ 
fected by the intrinsic accuracy of EKT, the performance 
of GE2, and the choice of basis set. We will not discuss 
the accuracy of EKT in this work, but will instead focus 
on the latter two points. To assess the performance of 
GE2 and GiE2, we have carried out a series of calcu¬ 
lations on several closed shell atoms and molecules. We 
start from atomic calculations since they are simpler, and 
then turn our discussion to small molecules. We have 
carried out UCCSD(T) calculations on each system to 
be used as a reference point throughout our discussion. 
We note that for the closed shell atoms and many of the 
molecules studied here the EA will be negative, meaning 
the system does not bind an extra electron, and in the 
complete basis set limit the EA would approach zero. 
However we include these systems in our analysis as a 
proof of concept, because we find the EA from GE2 with 
EKT for most cases agrees reasonably well with the re¬ 
sults from UCCSD(T) energy differences, as well as from 
HE Koopmans’ Theorem. 


A. Atoms 

In Eigure we present for a series of closed-shell 
atoms the IPs and EAs calculated with GE2 and GiE2 
using EKT, with HE using standard Koopmans’ Theo¬ 
rem, and with UCCSD(T) using energy differences, com¬ 
pared against experimental values when available. Eig¬ 
ure |2] illustrates that the IPs obtained from GE2 for these 
atoms are well converged with the basis set and sys¬ 
tematically underestimate experimental and UCCSD(T) 
IPs. In comparison with experimental values, the IPs 
can differ by up to 1 eV. Eor these atoms the best agree¬ 
ment occurs in the case of He, with a difference from 
experiment of around 0.1 eV. Eor comparison, calcula¬ 
tions were carried out non self-consistently (GiE2), and 
it was found that these values were in closer agreement 
with both UCCSD(T) and experimental IPs than the self- 
consistent GE2 IPs. This effect was observed previously 
in the work of Dahlen and van Leenwenp!8l [T9] . In con¬ 
trast, self-consistency appears to have a small effect on 
the EAs of these systems. The largest difference in elec- 
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GF 2 

GiF2 

HF 

UGGSD(T) 

Expt 

Be 2 

6.21 

6.99 

6.62 

7.42 


BHa 

12.82 

13.13 

13.52 

13.17 


C 2 H 2 

10.24 

11.24 

11.19 

11.36 

11.41 

C 2 H 4 

9.54 

10.22 

10.21 

10.55 

10.51 

CO 

12.20 

14.46 

15.09 

13.80 

14.01 

CO 2 

11.71 

12.88 

14.82 

13.61 

13.78 

H 2 CO 

9.12 

9.74 

12.02 

10.74 

10.88 

H 2 O 

11.31 

11.47 

13.86 

12.54 

12.65 

H 2 O 2 

9.51 

10.32 

13.31 

11.46 

11.70 

HGN 

12.26 

13.49 

13.50 

13.62 

13.61 

HF 

14.68 

14.66 

17.68 

16.02 

16.06 

Li 2 

4.69 

5.28 

4.95 

5.23 


LiF 

9.89 

9.70 

12.91 

11.37 


LiH 

7.77 

7.91 

8.20 

7.94 

7.90 

MgH 2 

9.80 

9.93 

10.09 

9.76 


N 2 

13.53 

14.97 

17.26 

15.36 

15.58 

Na 2 

4.67 

4.80 

4.51 

4.85 


NaF 

8.38 

8.15 

11.59 

9.98 


NaH 

6.82 

7.08 

7.43 

7.04 


NaLi 

4.69 

4.96 

4.71 

5.01 


NaOH 

6.42 

6.35 

9.11 

7.86 


NHs 

9.86 

10.10 

11.67 

10.76 

10.82 


TABLE L Calculated ionization potentials in eV using self-consistent GF2 and non self-consistent GiF2. Shown for com¬ 
parison are the ionization potentials from Hartree-Fock Koopmans’ Theorem and UCCSD(T) energy-differences, along with 
experimental values. All calculations use aug-cc-pVDZ. 


tron affinity occurs for Ne, with a difference of 0.2 eV. 
The majority of the atoms have a difference on the order 
of 1 meV between self-consistent and non self-consistent 
calculations. It should be noted that the EAs are not 
converged with respect to the basis set. Between basis 
sets the EAs can vary by around 1 eV, in a similar fashion 
to the UCCSD(T) electron affinities. 


B. Molecules 

Turning now to the molecules, in Tables [l| and [ll] we 
present our calculated IPs and EAs. Focusing first on 
the IPs, similar to the atoms, we find that GF 2 tends to 
systematically underestimate experimental values. Fur¬ 
thermore, comparing against UCCSD(T) reference val¬ 
ues in Fig. we find that GF 2 systematically under¬ 
estimates UCGSD(T) as well. In contrast, non-self- 


consistent GiF 2 again yields overall slightly larger IPs 
that are in better agreement with UGCSD(T) and exper¬ 
imental values than self-consistent GF 2 . This same effect 
has been observed in the work of Rostgaard, Jacobsen, 
and Thvgesen[28] with self-consistent GW and non-self- 
consistent GqWo applied to a similar set of molecular 
systems, where it was interpreted as being caused by 
overscreening in the former case and underscreening in 
the latter. GF 2 does not have a series of bubble dia¬ 
grams that are commonly thought to be responsible for 
the screening effects. However, the self-consistent GF 2 
has a series of composite type diagrams resulting from 
the iterative procedure where four of the diagrams from 
Figure are joined together into series of ladder like di¬ 
agrams. These diagrams most likely have effects similar 
to the series of bubbles in self-consistent GW which are 
causing over screening. We believe that since the self- 
consistent GF 2 redefines the Fock matrix during itera- 





FIG. 3. HF, GiF2, and GF2 IPs for molecules compared 
against reference UGGSD(T) values, calculated with aug-cc- 
pVDZ. 


tions, the results of the iterative procedure are less likely 
to overestimate the amount of correlation as in the case 
of MP2. Therefore, we assume that the good agreement 
of GiF 2 with UCCSD(T) is fortuitous. One could spec¬ 
ulate that in iterative GF 2 introducing the third order 
diagrams may be much more beneficial and would lead 
to systematically convergent IP results. 

o 

Examining the EAs now, we find GE 2 tends to give 
results which are slightly lower than both HE and 
UCCSD(T) values. Eor a few of the molecules with posi¬ 
tive EAs, GE 2 and GiE 2 do not recover the proper sign. 
However we emphasize that the EA is a notoriously dif¬ 
ficult property to calculate, and from comparing our re¬ 
sults to those from EKT-MPn in Table 2 of Ref. [45] it is 
not uncommon for a method to occasionally predict the 
incorrect sign for even small closed-shell molecules. Inter¬ 
estingly, whether or not the calculations are carried out 
self-consistently does not appear to have such a drastic 
effect on the EAs, as the GE 2 , GiE 2 , and HE EAs tend 
to be quite similar to the UCCSD(T) values for many 
of the molecules. We think the simplest explanation for 
this is the following: in the imaginary time Green’s func¬ 
tion EKT approach we are using, the EAs are essentially 
determined by the “hole”-part of G(r) near r = 0+, and 
the IPs likewise by the “particle”-part near r = (3. Be¬ 
cause these systems are small and weakly correlated, it 
may be that the “hole” or virtual orbital space does not 
relax as much between HE, GiE 2 , and GE 2 , as does the 
“particle” or occupied orbital space. Regardless, we find 
it encouraging that EAs can systematically be recovered 
from simply the imaginary time Green’s function of the 
neutral system, without need for considering molecular 
anions. 



GF 2 

GiF2 

HE 

UCCSD(T) 

Be 2 

-0.30 

-0.28 

- 0.21 

0.34 

BHg 

-0.92 

-0.94 

- 0.88 

-0.15 

C 2 H 2 

- 1.02 

- 1.02 

- 1.02 

-0.99 

C 2 H 4 

-1.09 

- 1.10 

- 1.10 

- 1.22 

CO 

-2.34 

-2.33 

-2.15 

-1.81 

CO 2 

-1.46 

-1.47 

-1.49 

-2.23 

H 2 CO 

-0.94 

-0.90 

-0.90 

-1.19 

H 2 O 

-0.93 

-0.96 

-0.96 

-0.75 

H 2 O 2 

-1.05 

-1.08 

-1.08 

-1.16 

HCN 

-0.83 

-0.79 

-0.79 

-0.69 

HF 

-0.95 

-0.97 

-0.97 

-0.80 

Li 2 

- 0.10 

- 0.10 

-0.08 

0.32 

LiF 

0.27 

0.29 

0.29 

0.34 

LiH 

0.20 

0.20 

0.21 

0.30 

MgH 2 

-0.45 

-0.44 

-0.44 

-0.50 

N 2 

-2.81 

-4.02 

-3.39 

-2.60 

Na 2 

-0.05 

-0.05 

- 0.01 

0.36 

NaF 

0.43 

0.43 

0.44 

0.48 

NaH 

0.22 

0.23 

0.24 

0.32 

NaLi 

-0.07 

-0.07 

-0.04 

0.35 

NaOH 

0.32 

0.31 

0.31 

0.39 

NH 3 

-0.94 

-0.97 

-0.98 

-0.75 

TABLE II. 

Galculated 

electron 

affinities 

in eV using self- 


consistent GF2 and non self-consistent GiF2. Shown for com¬ 
parison are the electron affinities from Hartree-Fock Koop- 
mans’ Theorem, and UGGSD(T) energy-differences. All cal¬ 
culations use aug-cc-pVDZ. 


VI. CONCLUSIONS 

In this work we have investigated how reliably IPs and 
EAs can be calculated from the Extended Koopmans’ 
Theorem (EKT) with an imaginary time Green’s func¬ 
tion in a second order approximation (GE 2 ). In con¬ 
trast to prior EKT works that determined the EA indi¬ 
rectly as the IP of the anion [45]. in this work we calcu¬ 
lated both IPs and EAs directly from the imaginary time 
Green’s function of the neutral system alone. Overall, 
we find that self-consistent GE2 with EKT recovers IPs 
and EAs that are systematically smaller than UCCSD(T) 
energy-differences and experiment reference values. In¬ 
terestingly, non-self-consistent GiE 2 on a Hartree-Eock 
reference consistently gives slightly larger IPs and EAs 
than self-consistent GE 2 , similar to what has been found 
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with GW and GqWo for IPs [28]. Because GF2 is defined 
in terms of the bare Goulomb interaction rather than 
the screened interaction as in GW, this suggests that the 
cause of the systematic underestimation of quasiparticle 
spectra by self-consistent vs non-self-consistent Green’s 
function methods may be more general than being specif¬ 
ically the result of over or underscreening caused by a se¬ 
ries of bubble diagrams present in the the GW approach. 

Regardless of the particular performance of GF2 or 
GiF2, the more general point of this work is that the 
EKT in conjunction with self-consistent Green’s function 
theory offers a reliable procedure for the unbiased theo¬ 
retical determination of quasiparticle spectra. Essentially 
the underlying scheme in this Green’s function EKT ap¬ 
proach is that the IPs are determined by the eigenvalues 
of the time-derivative of the “particle”-part of G(r) at 
T = /d, while the EAs are likewise found from the “hole”- 
part of G(r) at r = 0+. In this way the full quasi¬ 


particle spectra can in principle be reconstructed from 
simply the Green’s function on the imaginary time do¬ 
main, without the need for analytic continuation [3Qll32] 
or other numerical methods. Therefore the EKT allows 
one to obtain real axis quasiparticle spectra while still en¬ 
joying the computational benefits of using an imaginary 
axis Green’s function implementation. Eurthermore, the 
advantages of the EKT are that it is simple, is applicable 
to a Green’s function from any level of theory, requires 
only a trivial amount of computational time, and can be 
implemented in a blackbox manner. 
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